module dg_tools
implicit none
    ! 调整为更高精度（至少20位有效数字，指数范围至少300）
    integer, parameter :: qp = selected_real_kind(20, 300)
    
    ! 高精度常量（25位有效数字）
    real(kind=qp), parameter, dimension(4),private :: gauss_points = [ &
        0.3399810435848562648026658_qp, &
        0.8611363115940525752239465_qp, &
        -0.3399810435848562648026658_qp, &
        -0.8611363115940525752239465_qp ]
    
    real(kind=qp), parameter, dimension(4),private :: gauss_coef = [ &
        0.6521451548625461426269361_qp, &
        0.3478548451374538573730639_qp, &  
        0.6521451548625461426269361_qp, &
        0.3478548451374538573730639_qp ]   
    
    ! 使用高精度π值（至少25位）
    real(kind=qp), parameter :: pi = 3.141592653589793238462643_qp
contains
    ! 输出子程序（保持不变）
    subroutine show_consts()
        print '(A)', "Gauss_coefficient:"
        print '(F25.20)', gauss_coef  ! 调整输出格式为25位
        print '(A)', "Gauss_points :"
        print '(F25.20)', gauss_points ! 调整输出格式为25位
    end subroutine show_consts
    
    ! 获取子程序（保持不变）
    subroutine get_points(points)
        implicit none
        real(kind=qp), intent(out), dimension(4) :: points
        points = gauss_points
    end subroutine get_points
    subroutine get_pi(ppi)
        real(kind=qp), intent(out) :: ppi
        ppi = pi
    end subroutine get_pi
    ! 计算基函数在高斯点处的值（保持不变）
    subroutine compute_bvalue(k, b_value)
        implicit none
        integer, intent(in) :: k
        real(kind=qp), intent(out), dimension(k+1,4) :: b_value
        integer :: n
        select case (k)
            case (1)
                do n = 1,4
                    b_value(1,n) = 1
                    b_value(2,n) = gauss_points(n)
                end do
            case (2)
                do n = 1,4
                    b_value(1,n) = 1
                    b_value(2,n) = gauss_points(n)
                    b_value(3,n) = 0.5*(3*(gauss_points(n)**2) - 1)
                end do
        end select
    end subroutine compute_bvalue
    subroutine compute_dbvalue(k, db_value)
        implicit none
        integer, intent(in) :: k
        real(kind=qp), intent(out), dimension(k+1,4) :: db_value
        integer :: n
        select case (k)
            case (1)
                do n = 1,4
                    db_value(1,n) = 0
                    db_value(2,n) = 1
                end do
            case (2)
                do n = 1,4
                    db_value(1,n) = 0
                    db_value(2,n) = 1
                    db_value(3,n) = 3*gauss_points(n)
                end do
        end select
    end subroutine compute_dbvalue
    ! 积分计算子程序（保持不变）
    subroutine compute_int(fun_value,int_value)
        real(kind=qp), intent(in), dimension(4) :: fun_value
        real(kind=qp), intent(out) :: int_value
        int_value = dot_product(gauss_coef, fun_value)
    end subroutine compute_int

    !基函数边界点值获取
    subroutine compute_basis(k,basis)
        integer,intent(in) :: k
        real(kind=qp),intent(out), dimension(k+1,3) :: basis
        basis(1,:) = 1.0_qp
        basis(2,1) = -1.0_qp
        basis(2,2) = 1.0_qp
        basis(2,3) = 0.0_qp
        if(k==2) then
            basis(3,:) = 1.0_qp
        end if
    end subroutine compute_basis
    subroutine compute_dbasis(k,dbasis)
        integer,intent(in) :: k
        real(kind=qp),intent(out), dimension(k+1,3) :: dbasis
        dbasis(1,:) = 0.0_qp
        dbasis(2,:) = 1.0_qp
        if(k==2) then
            dbasis(3,1) = -3.0_qp
            dbasis(3,2) = 3.0_qp
            dbasis(3,3) = 0.0
        end if
    end subroutine compute_dbasis
    !生成质量矩阵
    subroutine mass_matrix(k,h,Mass)
        real(kind=qp), intent(in) :: h
        integer, intent(in) :: k
        real(kind=qp), intent(out) ,dimension(k+1,k+1) :: Mass
        real(kind=qp) :: c
        c = 2.0_qp / h
        Mass(1,1) = c*0.5_qp
        Mass(2,1) = 0.0_qp
        Mass(1,2) = 0.0_qp
        Mass(2,2) = c*1.5_qp
        if (k == 2) then
            Mass(1,3) = 0.0_qp
            Mass(2,3) = 0.0_qp
            Mass(3,3) = c*2.5_qp
            Mass(3,1) = 0.0_qp
            Mass(3,2) = 0.0_qp
        end if
    end subroutine mass_matrix
    ! 生成刚度矩阵
    subroutine stiff_matrix(k,stiff)
        integer, intent(in) :: k
        real(kind=qp), intent(out), dimension(k+1,k+1) :: stiff
        stiff(:,1) = 0.0_qp
        stiff(:,2) = 0.0_qp
        stiff(2,1) = 2.0_qp
        if (k==2) then
            stiff(:,3) = 0.0_qp
            stiff(3,2) = 2.0_qp
        end if
    end subroutine stiff_matrix
    !对初值做L2投影
    subroutine get_initial_coef(k,cell_num,h,u_value,a_coef,mass)
        integer,intent(in) :: k, cell_num
        real(kind=qp),intent(in) :: h
        real(kind=qp),intent(in), dimension(cell_num,4) :: u_value
        real(kind=qp),intent(out), dimension(k+1, cell_num) :: a_coef
        real(kind=qp),intent(in), dimension(k+1,k+1) :: mass
        real(kind=qp), dimension(k+1,4) :: b_value
        real(kind=qp), dimension(4) :: f_value
        real(kind=qp), dimension(k+1) :: right_b
        integer :: n,i
        call compute_bvalue(k, b_value)
        do n =1 , cell_num
            do i =1 ,k+1
                call vector_mul(4,b_value(i,:),u_value(n,:), f_value)
                call compute_int(f_value, right_b(i))
                right_b(i) = right_b(i) * (h/2.0_qp)
            end do
            a_coef(:,n) = matmul(mass, right_b)
        end do
    end subroutine get_initial_coef
    ! 向量对应元素相乘
    subroutine vector_mul(n,a,b,vec)
        integer, intent(in) :: n
        real(kind=qp), intent(in), dimension(n) :: a
        real(kind=qp) ,intent(in), dimension(n) :: b
        real(kind=qp), intent(out), dimension(n) ::vec
        integer :: i
        do i =1 ,n
            vec(i) = a(i) * b(i)
        end do
    end subroutine vector_mul
    ! 误差计算-L2范数
    subroutine compute_l2error(h,k,cell_num,a_coef,u_value,error)
        integer,intent(in) ::k, cell_num
        real(kind=qp),intent(in) :: h
        real(kind=qp),intent(out) :: error
        real(kind=qp),intent(in), dimension(k+1, cell_num) :: a_coef
        real(kind=qp),intent(in), dimension(cell_num,4) :: u_value
        real(kind=qp),dimension(k+1,4) :: b_value
        real(kind=qp),dimension(cell_num,k+1) :: a_trans
        real(kind=qp),dimension(cell_num,4) :: num_value
        real(kind=qp),dimension(4) :: tmp
        real(kind=qp),dimension(4) :: error_square
        real(kind=qp) :: tmp_int
        integer :: i
        call compute_bvalue(k, b_value)
        a_trans = transpose(a_coef)
        num_value = matmul(a_trans, b_value)
        do i=1, cell_num
            tmp = num_value(i,:) - u_value(i,:)
            call vec_square(4, tmp, error_square)
            call compute_int(error_square,tmp_int)
            tmp_int = (h/2.0_qp)*tmp_int
            error = error + tmp_int
        end do
        error = sqrt(error)
    end subroutine compute_l2error
    ! 计算误差-L1 范数
    subroutine compute_l1error(h,k,cell_num,a_coef,u_value,error)
        integer,intent(in) ::k, cell_num
        real(kind=qp),intent(in) :: h
        real(kind=qp),intent(out) :: error
        real(kind=qp),intent(in), dimension(k+1, cell_num) :: a_coef
        real(kind=qp),intent(in), dimension(cell_num,4) :: u_value
        real(kind=qp),dimension(k+1,4) :: b_value
        real(kind=qp),dimension(cell_num,k+1) :: a_trans
        real(kind=qp),dimension(cell_num,4) :: num_value
        real(kind=qp),dimension(4) :: tmp
        real(kind=qp) :: tmp_int
        integer :: i
        call compute_bvalue(k, b_value)
        a_trans = transpose(a_coef)
        num_value = matmul(a_trans, b_value)
        do i=1, cell_num
            tmp = abs(num_value(i,:) - u_value(i,:))
            call compute_int(tmp,tmp_int)
            tmp_int = (h/2.0_qp)*tmp_int
            error = error + tmp_int
        end do
        error = sqrt(error)
    end subroutine compute_l1error
    ! 计算阶数
    subroutine compute_order(n,error)
        integer, intent(in) :: n
        real(kind=qp), intent(in), dimension(n) :: error
        integer :: i
        real(kind=qp) :: order
        order = 0.0_qp
        do i = 1, n-1
            if (error(i+1) /= 0.0_qp) then
                order = log(error(i)/error(i+1)) /log(2.0_qp)
            end if
            print '(A,F12.8)', "Order of convergence: ", order
        end do
    end subroutine compute_order
    subroutine vec_square(n, vec, result)
        integer, intent(in) :: n
        real(kind=qp), intent(in), dimension(n) :: vec
        real(kind=qp), intent(out), dimension(n) :: result
        integer :: i
        do i = 1, n
            result(i) = vec(i)**2
        end do
    end subroutine vec_square
    subroutine basicCL_one_step_Euler(k, cell_num, a_coef, dt,basis,mass, stiff)
        integer,intent(in) :: k, cell_num
        integer :: nn, l, p
        real(kind=qp), intent(in) ::dt
        real(kind=qp), intent(in), dimension(k+1,k+1) :: mass, stiff
        real(kind=qp), intent(in), dimension(k+1,3) :: basis
        real(kind=qp), intent(inout), dimension(k+1,cell_num) :: a_coef
        real(kind=qp), dimension(cell_num,2) :: u_bound
        real(kind=qp), dimension(k+1) :: b
        !计算边界值
        u_bound = matmul(transpose(a_coef),basis(:,1:2))
        do nn = 1, cell_num
            if (nn ==1 ) then
                p = cell_num
            else
                p = nn - 1
            end if
            do l = 1, k+1
                b(l) = -u_bound(nn,2)*basis(l,2) + u_bound(p,2)*basis(l,1)
            end do
            a_coef(:,nn) = a_coef(:,nn) + dt * matmul(mass, b) + dt* matmul(mass, matmul(stiff, a_coef(:,nn)))
        end do
    end subroutine basicCL_one_step_Euler
    ! Rk3
    subroutine basicCL_RK3(k, cell_num, a_coef, dt, basis, mass, stiff)
        integer,intent(in) :: k, cell_num
        real(kind=qp), intent(in) :: dt
        real(kind=qp), intent(in), dimension(k+1,k+1) :: mass, stiff
        real(kind=qp), intent(in), dimension(k+1,3) :: basis
        real(kind=qp), intent(inout), dimension(k+1,cell_num) :: a_coef
        real(kind=qp), dimension(k+1,cell_num) :: a_coef1, a_coef2, recover
        recover = a_coef
        call basicCL_one_step_Euler(k, cell_num, a_coef, dt, basis, mass, stiff)
        a_coef1 = a_coef
        call basicCL_one_step_Euler(k, cell_num, a_coef1, dt, basis, mass, stiff)
        a_coef2 = 0.75*recover + 0.25*a_coef1
        call basicCL_one_step_Euler(k, cell_num, a_coef2, dt, basis, mass, stiff)
        a_coef = (1.0_qp/3.0_qp) * recover + (2.0_qp/3.0_qp) * a_coef2
    end subroutine basicCL_RK3
end module dg_tools